﻿Imports System.ComponentModel
''' <summary>
''' 
''' </summary>
''' <remarks></remarks>
Public Class LOGFFT
    Inherits TimeSeriesFrameTransformer

    Private _LogFrequenciesIndex As Integer()
    Private _LogBins As Integer = 64
    Private _LogBase As Double = Math.E
    Private _SampleRate As Integer = 22050
    Private _FrameWidth As Integer = 512
    Private _MaxFrequency As Double = 2855.0
    Private _MinFrequency As Double = 318.0
    Private _Phase As FFT.Phase = FFT.Phase.Signal
    Private _FFT As New FFT

    Public Sub New()

        MyBase.FrameStep = FrameStep._64
        MyBase.FrameWidth = FrameWidth._512

        _FFT.PhaseMode = _Phase
        GenerateLogFrequencies()

    End Sub

    Protected Overloads Overrides Function OnTransformFrame(frame As Double()) As Double()

        Dim complexSignal As Double() = New Double(2 * _FrameWidth - 1) {}

        For j As Integer = 0 To _FrameWidth - 1
            complexSignal(2 * j) = frame(j)
            complexSignal(2 * j + 1) = 0
        Next

        _FFT.Complex(complexSignal, True)

        Dim bins As Double() = ExtractLogBins(complexSignal)

        Return bins

    End Function

    Protected Overrides Sub OnSampleRateChanged(sampleRate As Integer)
        _SampleRate = sampleRate
        GenerateLogFrequencies()
    End Sub

    Protected Overrides Sub OnFrameWidthChanged(width As Integer)
        _FrameWidth = width
        GenerateLogFrequencies()
    End Sub

    ''' <summary>
    '''   Logarithmic spacing of a frequency in a linear domain
    ''' </summary>
    ''' <param name = "spectrum">Spectrum to space</param>
    ''' <returns>Logarithmically spaced signal</returns>
    Public Function ExtractLogBins(ByVal spectrum As Double()) As Double()
        Dim sumFreq As Double() = New Double(_LogBins - 1) {}
        For i As Integer = 0 To _LogBins - 1
            Dim lowBound As Integer = _LogFrequenciesIndex(i)
            Dim hiBound As Integer = _LogFrequenciesIndex(i + 1)
            'If lowBound = hiBound Then hiBound += 1
            For k As Integer = lowBound To hiBound '- 1
                Dim re As Double = spectrum(2 * k)
                Dim im As Double = spectrum(2 * k + 1)
                re /= _FrameWidth '/ 2    '//normalize img/re part
                im /= _FrameWidth '/ 2   '//doesn't introduce any change in final image (linear normalization)
                sumFreq(i) += Math.Sqrt(re * re + im * im)
            Next
            sumFreq(i) /= (hiBound - lowBound) + 1
        Next
        Return sumFreq
    End Function

    ''' <summary>
    '''        * An array of WDFT [0, 2048], contains a range of [0, 5512] frequency components.
    '''        * Only 1024 contain actual data. In order to find the Index, the fraction is found by dividing the frequency by max frequency
    '''        
    '''   Gets the index in the spectrum vector from according to the starting frequency specified as the parameter
    ''' </summary>
    ''' <param name = "freq">Frequency to be found in the spectrum vector [E.g. 300Hz]</param>
    ''' <param name = "sampleRate">Frequency rate at which the signal was processed [E.g. 5512Hz]</param>
    ''' <param name = "spectrumLength">Length of the spectrum [2048 elements generated by WDFT from which only 1024 are with the actual data]</param>
    ''' <returns>Index of the frequency in the spectrum array</returns>
    ''' <remarks>
    '''   The Bandwidth of the spectrum runs from 0 until SampleRate / 2 [E.g. 5512 / 2]
    '''   Important to remember:
    '''   N points in time domain correspond to N/2 + 1 points in frequency domain
    '''   E.g. 300 Hz applies to 112'th element in the array
    ''' </remarks>
    Private Shared Function FreqToIndex(ByVal freq As Double, ByVal sampleRate As Integer, ByVal spectrumLength As Integer) As Integer
        ' freq/((float) sampleRate/2)
        Dim fraction As Double = freq / (sampleRate / 2) '(freq + 0.5) / (sampleRate / 2)
        'N sampled points in time correspond to [0, N/2] frequency range 
        Dim i As Integer = CInt(Math.Round((spectrumLength / 2 + 1) * fraction))
        'DFT N points defines [N/2 + 1] frequency points
        Return i
    End Function

    Private Sub GenerateLogFrequencies()
        Dim logMin As Double = Math.Log(_MinFrequency, _LogBase)
        Dim logMax As Double = Math.Log(_MaxFrequency, _LogBase)
        Dim delta As Double = (logMax - logMin) / _LogBins

        _LogFrequenciesIndex = New Integer(_LogBins) {}

        Dim accDelta As Double = 0
        For i As Integer = 0 To _LogBins
            '32 octaves
            Dim freq As Double = Math.Pow(_LogBase, logMin + accDelta)
            accDelta += delta
            'Find the start index in array from which to start the summation
            _LogFrequenciesIndex(i) = FreqToIndex(freq, _SampleRate, _FrameWidth)
        Next

        'timederivwindow(Window, _SampleRate, 1)

    End Sub

    <Description("Number of bins."), DefaultValue(64)>
    Public Property LogBins As Integer
        Get
            Return _LogBins
        End Get
        Set(value As Integer)
            If _LogBins <> value Then
                _LogBins = value
                GenerateLogFrequencies()
            End If
        End Set
    End Property

    <Description("Maximum Frequency"), DefaultValue(2855.0)>
    Public Property FrequencyMax As Double
        Get
            Return _MaxFrequency
        End Get
        Set(value As Double)
            If value <> _MaxFrequency Then
                If value <= _MinFrequency Then Throw New ArgumentException("MaxFrequency cannot be less than or equal to MinFrequency.")
                If value >= _SampleRate Then Throw New ArgumentException("MaxFrequency cannot be greater than the SampleRate.")
                _MaxFrequency = value
                GenerateLogFrequencies()
            End If
        End Set
    End Property

    <Description("Minimum Frequency"), DefaultValue(318.0)>
    Public Property FrequencyMin As Double
        Get
            Return _MinFrequency
        End Get
        Set(value As Double)
            If value <> _MinFrequency Then
                If value >= _MaxFrequency Then Throw New ArgumentException("MinFrequency cannot be greater than or equal to MaxFrequency.")
                If value <= 0 Then Throw New ArgumentException("MinFrequency must be greater than 0.")
                _MinFrequency = value
                GenerateLogFrequencies()
            End If
        End Set
    End Property

    <DisplayName("Phase"), Description("FFT phase mode."), DefaultValue(FFT.Phase.Signal)>
    Public Property PhaseMode As FFT.Phase
        Get
            Return _Phase
        End Get
        Set(value As FFT.Phase)
            If value <> _Phase Then
                _Phase = value
                _FFT.PhaseMode = _Phase
            End If
        End Set
    End Property

    Public Overrides ReadOnly Property Description As String
        Get
            Return "Log Fast Fourier Transform"
        End Get
    End Property

End Class
